Replacement Class - Cleaning up and Spatial Joins

This notebooks will help me talk through a bit of material regarding spatial data, with some introduction to spatial joins.

Key topics for today:

• Using sf package
• Using tmap package
• Using tidycensus package
• Reminder on joins (focus: left join)
• Our Spatial Data
  • Neighborhoods
  • Joining with non-spatial data
  • Census data
  • Joining with spatial data
• ignore html in git (if time)

Packages

Standards:

library(knitr)
library(tidyverse)
library(janitor)
library(lubridate) # because we will probably see some dates
library(here) # a package I haven't taught you about before that doesn't do much, but ....
library(rnaturalearth)
library(WDI)
library(tigris)
library(rgdal)

Loading required package: sp
Please note that rgdal will be retired during 2023,
plan transition to sf/stars/terra functions using GDAL and PROJ
at your earliest convenience.
See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
rgdal: version: 1.6-2, (SVN revision 1183)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/gdal
GDAL binary built with GEOS: FALSE 
Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.5-1
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

library(sp)

Some additional packages focused on today’s work:

library(sf) # working with simple features - geospatial

Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(tmap)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(tidycensus)

Informational resources

• An overall resource on mapping in R: https://bookdown.org/nicohahn/making_maps_with_r5/docs/introduction.html
• A starting point to learn about sf: https://r-spatial.github.io/sf/articles/
• Getting started with tmap: https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html
• The tidycensus package: https://walker-data.com/tidycensus/index.html
• The book on tidycensus : https://walker-data.com/census-r/index.html

Using the Neighborhood Geospatial Data (using /data)

Our first data source comes from opendata.dc

https://opendata.dc.gov/datasets/DCGIS::dc-health-planning-neighborhoods/about

I will use the GeoJSON file. (Newer, not necessarily better, but … a single file. Not smaller, but … this one is not big.)

Data is easily readable

neigh=st_read(here("DC_Health_Planning_Neighborhoods_joey.geojson")) %>% clean_names()

Reading layer `DC_Health_Planning_Neighborhoods' from data source 
  `/Users/joewhite/Desktop/R-Studio/ds241/work/DC_Health_Planning_Neighborhoods_joey.geojson' 
  using driver `GeoJSON'
Simple feature collection with 51 features and 8 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.11976 ymin: 38.79165 xmax: -76.9094 ymax: 38.99556
Geodetic CRS:  WGS 84

class(neigh)

[1] "sf"         "data.frame"

plot(neigh)

Reminder - Joins

df1=tibble(fruit=c("apple","banana","cherry"),cost=c(1.5,1.2,2.25))
df2=tibble(fruit=c("apple","apple","cherry","lemon"),
           desert=c("pie","cobbler","cobbler","cheesecake"),
           cal=c(400,430,500,550))
df1

df2

left_join(df1,df2,by="fruit")

Investigating joining spatial and non-spatial data

Covid case information is available from opendatadc:

https://opendata.dc.gov/datasets/DCGIS::dc-covid-19-total-positive-cases-by-neighborhood/about

Read cases information:

df_c=read_csv(here("DC_COVID-19_Total_Positive_Cases_by_Neighborhood_joey.csv")) %>% clean_names() 

df_cases=df_c %>%
  filter(as_date(date_reported) == "2022-02-22") %>% 
  separate(neighborhood,into=c("code","name"),sep = ":") %>%
  mutate(code=case_when(
    code=="N35" ~"N0",
    TRUE ~ code
  )) %>%
  select(-objectid,-date_reported)

Regular joining (of dataframes)

neigh2=left_join(neigh,df_cases,by=c("code")) 

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(neigh2) +tm_polygons("total_positives",alpha=.5)

Joining with other spatial data

Let’s get some data using tidycensus. Need an API key https://api.census.gov/data/key_signup.html

census_api_key("d44395e2fa101f82260fae6b845676d71f017b70")

To install your API key for use in future sessions, run this function with `install = TRUE`.

#what variables
v20 = load_variables(2018,"acs5")
# median_family_income="    B06011_001" 
# all "B00001_001"  
#black "B02009_001"

Get some data:

df_cencus=get_acs(geography = "tract",
                  variables=c("median_inc"="B06011_001",
                              "pop"="B01001_001",
                              "pop_black"="B02009_001"),
                  state="DC",geometry=TRUE,year=2018) 

Getting data from the 2014-2018 5-year ACS
Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
Using FIPS code '11' for state 'DC'

  |                                                                                          
  |                                                                                    |   0%
  |                                                                                          
  |==============================                                                      |  35%
  |                                                                                          
  |==================================                                                  |  40%
  |                                                                                          
  |==============================================                                      |  55%
  |                                                                                          
  |=================================================                                   |  58%
  |                                                                                          
  |================================================================                    |  76%
  |                                                                                          
  |=====================================================================               |  82%
  |                                                                                          
  |==============================================================================      |  93%
  |                                                                                          
  |====================================================================================| 100%

class(df_cencus)

[1] "sf"         "data.frame"

plot(df_cencus)

It’s in long format. Let’s make it wide.

df_cens=df_cencus %>% select(-moe) %>% spread(variable,estimate) 

tm_shape(df_cens) +tm_polygons("median_inc",alpha=.5)

  tm_shape(neigh2) +tm_borders(col="blue",lwd=5,alpha=.2)+
  tm_shape(df_cens) +tm_borders(col="red",lwd=1,alpha=.3)

#<<<<<<< HEAD
#df_j=st_join(df_cens,neigh2)
#=======
#df_j=st_join(df_cens,neigh2,prepared=FALSE)
#>>>>>>> aaf01be5cf721819dd2df615aef7a1999bcec0c2

df_cens_adj=df_cens %>% st_transform(4326)

df_j=st_join(df_cens_adj,neigh2,largest=TRUE)

Warning: attribute variables are assumed to be spatially constant throughout all geometries

Other order?:

#<<<<<<< HEAD
##df_j_rev = st_join(neigh2,df_cens_adj,largest=TRUE)
#=======
#df_j_rev = st_join(neigh2,df_cens_adj,largest=TRUE)
#>>>>>>> aaf01be5cf721819dd2df615aef7a1999bcec0c2

Since we want the geometry for the NEIGHBORHOODS, we need a to work a little harder:

df1=df_j %>% select(median_inc,pop,pop_black,objectid) %>%
  group_by(objectid) %>%
  summarise(pop_n=sum(pop),
            pop_black_n=sum(pop_black), 
            adj_median_income=sum(pop*median_inc)/pop_n) 

plot(df1)

#df2=left_join(neigh2,df1)

df2=left_join(neigh2,df1 %>% st_set_geometry(NULL))

Joining, by = "objectid"

df2=df2 %>% mutate(black_perc=pop_black_n/pop_n, covid_rate=total_positives/pop_n)
tm_shape(df2)+tm_polygons(c("adj_median_income","covid_rate","black_perc"))

df2 %>% filter(objectid!=30) %>% tm_shape()+tm_polygons(c("adj_median_income","covid_rate","black_perc"),alpha=.4)

#find where people ride bikes (bikes started in each ‘district’)

bike_data <- read_csv("202209-capitalbikeshare-tripdata.csv") %>% clean_names()

Rows: 386085 Columns: 13── Column specification ──────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (5): ride_id, rideable_type, start_station_name, end_station_name, member_casual
dbl  (6): start_station_id, end_station_id, start_lat, start_lng, end_lat, end_lng
dttm (2): started_at, ended_at
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

ggplot(bike_data, mapping = aes(x = start_lng, y = start_lat)) +
  geom_point()

bike_data_sf <- bike_data %>%
  mutate_at(vars(start_lat, start_lng), as.numeric) %>%
  st_as_sf(
    coords = c("start_lat", "start_lng"), 
    agr = "constant",
    crs = "4326"
  ) %>%
  rename(point_location = 'geometry') %>%
  sample_n(1000)

st_set_crs(4326) %>% st_transform(4326)

, agr = “constant”, crs = “4326” )

st_crs(bike_data_sf)

Coordinate Reference System: NA

st_crs(neigh2)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

st_crs(bike_data_sf) <- 4326

df_plz = st_join(neigh2, bike_data_sf,largest=TRUE)